Cross-correlations in scaling analyses of phase transitions 
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Thermal or finite-size scaling analyses of importance sampling Monte Carlo time series in the 
vicinity of phase transition points often combine different estimates for the same quantity, such as 
a critical exponent, with the intent to reduce statistical fluctuations. We point out that the origin 
of such estimates in the same time series results in often pronounced cross-correlations which are 
usually ignored even in high-precision studies, generically leading to significant underestimation of 
statistical fluctuations. We suggest to use a simple extension of the conventional analysis taking 
correlation effects into account, which leads to improved estimators with often substantially reduced 
statistical fluctuations at almost no extra cost in terms of computation time. 
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With the advent of the renormalization group in the 
1960s, the notions of scaling and universality have been 
combined into the solid basis of our understanding of crit- 
ical phenomena in statistical physics, field theory and a 
wealth of applications in areas ranging from solid-state 
physics [H to cosmology Q . It is through the remarkable 
fact that the most important properties of a continuous 
phase transition are independent of many microscopic 
details and, instead, only depend on a small number of 
fundamental characteristics of a system, such as the di- 
mensionality and the symmetries of the order parameter, 
that we can accurately describe such different situations 
as, e.g., the liquid- vapor phase transition and the ferro- 
magnetic transition of uniaxial magnets, with one and the 
same scaling theory. Recently, the investigation of quan- 
tum phase transitions has opened up a new Pandora's 
box with a wealth of transitions partially defining novel 
universality classes [3] • The direct applicability of results 
from simple models to a range of experimentally realized 
systems implied by the principle of universality renders 
high-precision determinations of critical parameters for 
the most common universality classes a rewarding goal. 

Particularly through the conception of advanced finite- 
size scaling (FSS) approaches and novel efficient algo- 
rithms [4J, Monte Carlo (MC) simulations have grown 
up to become a tool for the determination of univer- 
sal critical quantities clearly competitive compared to 
the more traditional approaches of high-temperature and 
field-theoretic expansions Likewise, the detailed in- 
vestigation of systems undergoing first-order phase tran- 
sitions has become a classic application of the MC simu- 
lation technique Q. Major advances in the competitive- 
ness of the MC method for these purposes came with 
the advent of histogram methods [7( and generalized- 
ensemble simulation techniques such as the multicanon- 
ical method H, both of which allow for extracting es- 



timates of thermal averages for a continuous range of 
temperatures or other external parameters from a single 
MC simulation. It is only through this effective con- 
tinuity of information that high-precision investigations 
of phase transitions have come into the reach of simula- 
tion methods. For arriving at high-precision estimates, 
however, all possible sources of error must be put under 
close scrutiny. This is often done to a high degree of 
sophistication concerning the systematic errors resulting 
from corrections to scaling [B| and the statistical errors 
resulting from the stochastic nature of MC time series 
(including their timewise autocorrelations for the case of 
the most commonly used Markov chain MC techniques) 
[tj [l(j ■ It has not been systematically discussed previ- 
ously, however, that the extraction of different estimates 
from a single time series in thermal or FSS analyses must 
entail cross-correlations. As will be shown below, neglect- 
ing their effect not only results in systematically wrong 
estimates of statistical errors, but also fails to fully ex- 
ploit the available time-series data to yield the maximum 
statistical precision obtainable. 

Although our considerations apply generally to all sit- 
uations where a number of different estimates from the 
same (set of) simulation(s) are combined to a final re- 
sult, for specificity we consider the FSS analysis of sim- 
ulation data in the vicinity of a critical point. For the 
purpose of illustration we choose the technique outlined 
in Ref . [ll| , but very similar considerations apply to alter- 
native approaches, see, e.g., Refs. 0E3]. To be specific, 
we here use a magnetic language and first consider the 
maxima of the derivative of the magnetization cumulants 
U 2k = l- (\m\ 2k )/3(\m\ k ) 2 for k = 1, 2, . . .: 
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where L denotes the linear size of the system and (3 is the 
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inverse temperature. This relation allows for a precise de- 
termination of the correlation length exponent v without 
previous knowledge of the critical temperature. In many 
cases at least some of the scaling corrections, such as an 
effective leading correction with exponent w as indicated 
in Eq. ([T]), need to be taken into account to achieve the 
desired level of accuracy. An analogous relation holds for 
the scaling of the maxima of the logarithmic derivative 
of magnetization moments, 

dln(H fc ) 
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These scaling relations for determining v only become 
useful as soon as the maximum values (d£/ 2 fc/d/3)| max 
and (dln(|m| fe )/d/3)| max can be computed to high accu- 
racy without the need for repeated simulations manually 
tracking their locations in j3. In case of a histogram or 
reweighting analysis of a single canonical simulation, this 
is effected through the continuous family of estimates 
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for the thermal average (A)p from a time series {At} 
of measurements resulting from an importance sampling 
simulation at inverse temperature /3o- Conventional tech- 
niques of numerical analysis such as a golden section 
search then allow for an efficient determination of the 
maxima of Eqs. (JTJ) and Q to high precision. Once v 
has been determined, the scaling of the shifts of the lo- 
cation of the maxima of quantities such as d[/2fc/d/3 and 
dln(|m| fc )/d/3 as well as the specific heat, susceptibility 
etc. allow to locate the transition coupling (3 C . Finally, 
the remaining critical exponents may be estimated from 
the FSS of the maxima of the specific heat to yield a/v, 
of the susceptibility to yield 7/1/ etc. Since the exponent 
v enters all FSS relations, it clearly is of utmost impor- 
tance to exploit the available data to their fullest for a 
precise estimate of v. In view of the family of relations 
([T|) and @, this certainly includes a combination of es- 
timates from dt/ 2fc /d/3 and dln(|m| fe )/d/3 as well as from 
the different choices of the parameter k 

To see how this combination should be performed, con- 
sider a number n of different estimators X4 with the 
same expectation value (e.g., Xi = Pi). A combined 
average results from a linear combination x = J2i a i%i 
with Ei a « = 1- While any such combination yields a 
valid estimator of (x), e.g., the arithmetic mean x p i a i n 
with at = l/n, the ensuing statistical fluctuations will 
be larger than necessary. For uncorrelated estimates x$ 
minimal variance of x is achieved for the error-weighted 
mean x orr with [14 1 

where a 2 (xi) denotes the variance of Xj and Z = 
1/ <r 2 {xi). In general, however, the estimates Xj, stem- 



ming from a reweighting analysis of the same MC time 
series, will be substantially correlated. Under these cir- 
cumstances, the optimum choice is a covariance- weighted 
mean x cov with weights [3, EH 
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where T{x) 1 denotes the inverse of the covariance ma- 
trix T lj (x) = (xiXj) - (xi){Xj) and Z = EyPX^) -1 ]^'- 
Since for uncorrelated estimates [r(.T) _1 ]y = / '<j 2 {xi), 
Eq. ^ is recovered in this special case. Even more dra- 
matically affected by correlations are the statistical er- 
rors of averages, where the standard formula Cuncorr(cc) = 
~^2iOi 2 (J 2 {xi) is no longer valid and must be modified to 
read a 2 orr (x) = J^i j cti&jFij (x), generically leading to an 
underestimate of fluctuations via the naive (and wrong) 
estimator cr 2 ncm . r . 

To check for the strength of such correlation effects and 
their influence on finding optimal averages endowed with 
valid estimates of statistical errors, we performed a FSS 
analysis of the critical points of the ferromagnetic Ising 
model in two (2D) and three (3D) dimensions. Time 
series data for the configurational energy and magnetiza- 
tion were produced from one single-cluster update sim- 
ulation [4| per system size at a fixed temperature. Esti- 
mates for the exponent v were extracted from FSS fits of 
the relations |1} resp. to the maxima of d£/2fe/d/3 with 
k = 1 and k = 2 as well as dln(|m| fc )/d/3 with k = 1, 2, 
and 3 extracted from a reweighting analysis. Statistical 
errors for the individual estimates were calculated via 
a jackknife analysis Q over the reweighting procedure, 
taking timewise autocorrelations into account. Likewise, 
the covariance matrix T was determined from the non- 
parametric jackknife estimator known to be especially 
robust @, 
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Here, n denotes the number of jackknife blocks, Vi( s ) de- 
notes the value for jackknife block s and is the arith- 
metic average of the Vi( s )- F° r the results presented here, 
n = 100 blocks were used, where we checked that the re- 
sults are invariant, at the level of statistical fluctuations, 
to the choice of a significantly larger number of blocks. 

For the case of the 2D model, simulations were per- 
formed at the asymptotic critical coupling j3 c — h ln(l + 
V2) = 0.440686794 . . ., using a series of square lattices 
of linear size L = 16, 24, . . ., 192. For this model, 
and the considered range of system sizes, we do not find 
corrections to scaling to be very pronounced, such that 
high-quality fits can be achieved while ignoring the terms 
proportional to L~ w in (fTJ and ([2]) and restricting the 
range of system sizes to L > 32. The resulting estimates 
are collected in the left part of Table |TJ Table fll] shows 
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TABLE I: Estimates of the correlation length exponent v for 
the 2D and 3D Ising models from the scaling of the maxima 
(0 and ((2|, as well as different averages and error estimates 
explained in the main text. The 2D reference value is exact, 
whereas in 3D it is taken from the recent review [f|. 





2D 




3D 






V 


a 


v 


a 


dln(H) 

dB 
dln{m 2 ) 

dB 
dhi(|m.| 3 ) 


1.0085 
1.0128 


0.0183 
0.0194 


0.6358 
0.6340 


0.0127 
0.0086 


1.0175 


0.0201 


0.6326 


0.0062 


dU 2 


1.0098 


0.0281 


0.6313 


0.0020 


dS 4 
dp 


1.0149 


0.0511 


0.6330 


0.0024 


X plain ^"uncorr 


1.0127 


0.0141 


0.6334 


0.0038 


V corr 




0.0269 




0.0067 


X CTT O uncorr 


1.0123 


0.0102 


0.6322 


0.0015 


& corr 




0.0208 




0.0024 


^cov 0" co rr 


0.9935 


0.0078 


0.6300 


0.0017 


reference value 


1 




0.6301 


0.0004 



the matrix of correlation coefficients pij = Tij/uiUj for 
these estimates as computed from the jackknife approach 
([6]) . ft maybe does not come unexpected that all of the 
estimates for zy, resulting from structurally similar ob- 
servables in the magnetic sector, are highly correlated 
with pij > 0.8. One might naturally wonder, then, if 
it is indeed worthwhile to consider all of these differ- 
ent estimators instead of, say, the single most precise 
one. The different averages discussed above are listed 
in the lower part of Table Q] together with the error es- 
timates cr uncorr neglecting correlations and <t C oit taking 
them into account. For the plain average as well as 
the error- weighted mean it is apparent that, although 
funcorr seems to indicate smaller fluctuations than for 
any single estimate, using the proper error er corr the sit- 
uation is reversed and the uncertainties of some of the 
single estimates, namely those stemming from the loga- 
rithmic derivatives, are smaller than the true fluctuation 
of these averages. For the full covariance- weighted mean, 
on the other hand, one arrives at v cov = 0.9935(78), 
which has clearly smaller fluctuations than any of the 
individual estimates. As is apparent from the lower part 
of Table [Til this improvement is effected through a dra- 
matically different choice of weights for the individual 
estimates as compared to the error-weighting or plain- 
average schemes. Comparing the standard deviations of 
the most commonly used average v elT and the new z> cov , 
it is striking that statistical precision is increased by al- 
most a factor of three merely by using different weights in 
the average. Against our usual intuition developed from 
statistics of uncorrelated events, the average v cov is here 
found to be smaller than all individual estimates. This 



TABLE II: Correlation coefficients pij = Vij joi(Jj between 
estimates of the critical exponent v of the 2D Ising model 
extracted from the maxima JT} and ((2]). The lower part of 
the table shows the weights at of the individual estimates in 
the plain, error-weighted and covariance-weighted averages, 
respectively. 



dln(|mj) dln(m 2 } dln(|m|' 



dln(|m|) 

df3 
dln(m 2 > 

dj3 
dln(|m| 3 ) 

dB 
dth 

dp 
dU 4 

dp 



dp 



dp 



dp 



dU 2 
dp 



dUj 
dp 



1.000 
0.974 
0.939 
0.920 
0.897 



0.974 
1.000 
0.991 
0.817 
0.869 



0.939 
0.991 
1.000 
0.743 
0.820 



0.920 0.897 

0.817 0.869 

0.743 0.820 

1.000 0.860 

0.860 1.000 



<^z, plain 


1.000 


1.000 


1.000 


1.000 


1.000 




0.315 


0.271 


0.248 


0.034 


0.132 


Q-i ,cov 


5.007 


-2.426 


-0.281 


-0.104 


-1.196 



is illustrated in Fig. [TJ where also be interpreted 

as a correlated fit to a constant (see also Ref. fl6j|). 

Simulations of the 3D Ising model were performed for 
simple cubic lattices of edge lengths L = 8, 12, 16, 
. . ., 128 at the fixed coupling f3 = 0.22165459 reported 
in a high-precision study as estimate for the transition 



point [17j . Here, scaling corrections for the logarith- 
mic derivatives of magnetization moments are sufficiently 
pronounced to warrant the inclusion of the L~ w correc- 
tion term of Eq. ((2]). For the cumulants, corrections are 
so small that, instead, fits of the uncorrected form were 
used on the range L > 32. The resulting estimates of v 
are collected on the right side of Table UJ Concerning the 
various averages, it is again found that errors are clearly 
underestimated when neglecting correlations, and for the 
plain and error-weighted means the true fluctuations are 
indeed larger than the errors of the single most precise es- 
timate. In contrast, the covariance-weighted mean yields 
^cov = 0.6300(17), significantly more precise than the 
single estimates as well as the averages not taking corre- 
lations into account. 

Similar considerations apply to the correlations be- 
tween estimates of different exponents. In particular, 
taking the scaling relations for critical exponents into 
account, the magnetic and energetic scaling dimensions 
might be estimated from different observables. For in- 
stance, the magnetic scaling dimension can be estimated 
via the relations x a = (3/u and x a = d/2 — ^/2v from the 
FSS of the magnetization at its inflection point and the 
magnetic susceptibility at its maximum via the relations 
(|m|) inf (L) = m -L _/3/l/ and Xmax(i) = XoL l/l/ , respec- 
tively. Table IIIII summarizes the correlation analysis for 
x a in the 2D model, where through the pronounced anti- 
correlation of the two estimates of x a the uncorrelated 
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FIG. 1: (Color online) Estimates ol the critical exponent v of 
the 2D Ising model from the FSS of the indicated observables. 
The horizontal lines indicate the plain average v v \^ 1TL (dotted), 
the error- weighted mean v e „ (dashed), and the covariance- 
weighted mean j? cov (solid). The shaded areas indicate the 
corresponding one-sigma environments v ± cr C orr • 



error er U ncorr over-estimates statistical fluctuations, and 
already the error-weighted mean is somewhat more pre- 
cise than either of the two single estimates. Still, the 
covariance-weighted mean is even more precise, yielding 
S(T,cov = 0.125(1), directly at the exact value x a = 1/8. 
For the 3D model, on the other hand, (only) the correla- 
tion analysis reveals that both estimates of x a are nearly 
uncorrelated such that, for this specific case, the full re- 
sult approximately coincides with the naive approach ne- 
glecting correlations. 

To summarize, we have seen that substantial cross- 
correlations exist between quantities estimated via his- 
togram analyses from time series of Markov chain MC 
simulations. Taking these into account by a straightfor- 
ward extension of the common data analysis reveals a 
generic underestimation of statistical error by the con- 
ventional approach. On the other hand, it suggests im- 
proved estimators with often substantially reduced statis- 
tical fluctuation resulting, for some examples, in a three- 
fold reduction of statistical error which could otherwise 
only be achieved with an about tenfold increase of sim- 
ulation time with the conventional analysis. While these 
effects have been illustrated here for the case of the criti- 
cal exponents of the Ising model, very similar behavior is 
expected for non-universal quantities, including proper- 
ties of first-order transitions [6j , and for different appli- 
cations, including the problems in soft-matter systems 
[HI], for quantum critical points or the extremely 
costly simulations of disordered systems [l9j. These ap- 
plications, together with the flexibility in choosing dif- 
ferent thermal or FSS approaches, render the outlined 
technique quite generic. 

M.W. acknowledges support by the DFG through 
the Emmy Noether Programme under contract No. 
WE4425/1-1. 



TABLE III: Correlation analysis and averages of estimates 
of the magnetic scaling dimension x a of the 2D Ising model 
from the scaling of the magnetization at its inflection point 
and the magnetic susceptibility at its maximum. 



fits 



corr. coeff. /weights 







a 


(H)inf 


Xmax 


(M)inf 


0.1167 


0.0054 


1.0000 


-0.6414 


Xmax 


0.1271 


0.0020 


-0.6414 


1.0000 


^plain ^uncorr 


0.1219 


0.0027 


1.0000 


1.0000 


<^corr 




0.0021 






^crr ^uncorr 


0.1261 


0.0016 


0.0944 


0.9056 


& corr 




0.0013 






%cov <^corr 


0.1250 


0.0010 


0.2050 


0.7950 


reference value 


0.125 
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